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Abstract 

In this work we study, at the single molecular level, the thermodynamic and dynamic character- 
istics of an enzymatic reaction comprising a rate limiting step. We investigate how the stability of 
the enzyme-state stationary probability distribution, the reaction velocity, and its efficiency of en- 
ergy conversion depend on the system parameters. We employ in this study a recently introduced 
formalism for performing a multiscale thermodynamic analysis in continuous-time discrete-state 
stochastic systems. 
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I. INTRODUCTION 



The vast majority of processes taking place inside a cell consist of, or at least involve 
to some extent, chemical reactions. Moreover, most, if not all of the chemical reactions 
taking place inside a cell are catalyzed by enzymes. Therefore, a profound understanding of 
enzymes' performance is necessary to better comprehend the processes of life. 

Homeostasis, understood as the coordinated physiological processes which maintain most 
of the steady states in an organism, is regarded as a landmark concept in biology [1] . It can be 
found, to some extent, in all living beings, and allows them to perform in optimal conditions 
despite ever-changing surroundings and inputs. From a dynamical standpoint, homeostasis 
implies the existence of a stable steady state [2H1]. Thus, the quality of homeostasis can 
be measured by the volume of the steady-state basin of attraction in phase space and/or 
the relaxation time with which the system returns to the steady state after a perturbation. 
Having a large basin of attraction is important because it allows the system to come back to 
the steady state even in the face of large deviations. On the other hand, a rapid relaxation 
time permits the system to quickly recover an optimal state after it is perturbed. 

Recent studies on finite-time thermodynamic engines and heat pumps have shown that 
their stability and their thermodynamic performance are often governed by the same param- 
eters. It has been observed that the system stability usually weakens as its thermodynamic 
properties improve [oTfPT]. Consequently, these parameters need to be tuned to achieve an 
optimal trade-off between favorable thermodynamic and dynamic properties. Similar studies 
on the stretch-reflex regulatory pathway and on a simple Brownian motor have confirmed 
these findings jTHJ [19] , in agreement with the notion that good design principles are usually 
shared by both artificial and biological systems [2]. If these results are of general applica- 
bility to a wide range of intracellular energy-converting processes, it would mean that the 
maintenance of the cell homeostatic state entails an expenditure of energy, which has to 
be taken into consideration to understand how organisms adapt to a constantly changing 
environment. 

In the present work we study, at the single molecular level, the thermodynamic and 
dynamic characteristics of an enzymatic reaction. For the sake of simplicity we advocate to 
enzymatic reactions comprising a rate limiting step. We investigate how the stability of the 
enzyme-state stationary probability distribution, the reaction velocity, and its efficiency of 
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energy conversion depend on the system parameters. We employ in this study a recently 
introduced formalism for performing a multiscale thermodynamic analysis in continuous- 
time discrete-state stochastic systems. 

II. MODELING ENZYMATIC REACTIONS 

A simple but comprehensive enough model for an enzymatic reaction consists in picturing 
the enzyme as undergoing the following series of chemical reactions [20]: first, a free enzyme 
E binds the substrate S; then, the bound substrate is converted into the product P to form 
the enzyme-product complex Ep\ and finally, the product is released leaving the enzyme 
free to catalyze another reaction. This process can be summarized as follows using the 
conventional notation for chemical reactions: 

E + S ^ E s ^ E P ^ E + P. 

Of all these reactions, the conversion of the substrate S into the product P (Es f± Ep) 
is in many cases the rate limiting process. Taking this into account and assuming that 
the substrate ([S]) and product ([P]) concentrations remain constant along the catalytic 
reaction, we can visualize a single enzyme as going through a series of transitions that 
change the enzyme state cyclically during the catalytic process; see Figure [1} There, the 
enzyme state is represented as where index i denotes the i th cycle — that in which the 
i th product molecule is synthesized — while j = 1,2,3 respectively correspond to states Ep, 
E + P, and E$. The assumption that [S] and [P] remain constant allows us to regard kj~ and 
kj (j = 1,2,3) as pseudo first order reaction rates. Finally, the assumption that E$ f± Ep 
is the rate limiting process implies that vertical transitions in the scheme of Figure [T] are 
much faster processes than those involving changes in index i. 

III. PROBABILISTIC DESCRIPTION AND TIME-SCALE SEPARATION 

Let us introduce a probabilistic description for an enzyme reaction at the single molecule 
level. And let us follow as well the approach in [211 122] to simplify the model, taking 
advantage of the assumed separation of time scales. Let P(i,j;t) denote the probability 
that the enzyme is in state at time t. From the scheme in Figure [ll the chemical 
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FIG. 1: Schematic representation of the various states available for an enzyme molecule and the 
transitions rates among them, while the enzyme is catalyzing the synthesis of molecules P. The 
states are denoted as where j = 1,2,3 respectively correspond to states Ep, E + P, and 

Es, while the index i denotes the i th enzyme cycle: that in which the i th product molecule is 
synthesized and released. 

master equation for P(i,j;t) consists of the following set of coupled differential equations: 
dP(i,l;t) 



dt 

dP(i,2;t) 
dt 

dP(i,3;t) 



k+P{i - 1, 3; t) + KP(i, 2; t) - + k+)P(i, 1; t), (1) 
k+P(t, 1; t) + k£P(i, 3; t) - + k+)P(i, 2; t), (2) 
k+P(i, 2; t) + k^P(t + 1, 1; t) - (k 2 + k+)P(z, 3; t). (3) 



dt 

The probability that the enzyme is in state (i, •) at time t is given by 

3 

P(z;t) = Y,P(hr,t). (4) 

3=1 

On the other hand, it follows from the definition of conditional probability that 

P{i,j-t) = P{j\i-t)P{i-t). (5) 

Add Equations ([[])- ([3]) and use Q and (JS]) to obtain 
dP(i; t) 



dt 



= (k+P(3\i-l;t))P(i-l;t) + (k 3 P(l\i + l;t))P(i + l;t) 
- {k+P(3\i;t) + k 3 _P(l\i;t)) P(i;t). (6) 



Differentiate ^ and assume a time-scale separation so that the transitions between states 
(i — 1, 1) and (i, 3) are much slower than all the other (that is, they are the rate limiting 
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dP(l\i 


t) 


dt 




dP(2\i 


t) 


dt 




dP(3\i 


t) 


dt 



steps along the reaction chain). Then, 

k^P{2\i;t) + kfP(l\i;t), (7) 

kfP(l\i; t) - k^P(2\i; t) + k^P(3\i; t) - fc^P(2|z; t), (8) 

k$P(2\i;t)-lqP(3\i;t). (9) 

If we invoke once more the separation of time scales to assume that the fast dynamics 
rapidly reach an equilibrium distribution while slow dynamics have not changed noticeably: 
dP(j\i; t)/dt m 0, j = 1, 2, 3, we have that 

P{l\i;t) ~ P*{l\i), 

P(2\i;t) ~ P*(2\i) = K x P*{l\i), 

P{2\i-t) ~ P*(3|z)=K 2 P*(2|z) = K 1 ^ 2 P*(l|i), 

where i^- = k^/kj, j = 1,2,3. Finally, the normalization condition (J^=i P*(j\i) = 1) 
implies that 

p *< 1|! > = uft'w < 10 > 

*w - mffk^- (12) 

Notice that P*(j\i) (j = 1,2,3) are all independent of %. Furthermore, it is straightforward 
to prove that this stationary distribution satisfies the following relations: 

k+P*(l\i) = k^P*(2\i), and k+P*(2\i) = k 2 P*(3\i). 

This last result implies that the fast dynamics are in equilibrium only if each one of the 



underlying chemical reactions is in equilibrium itself. Finally, substitution of (10)-(12) into 
(|6| allows us to conclude that, when time-scale separation is possible and the enzyme states 
are grouped as sketched in Figure [TJ the system dynamics is that of a biased one-dimensional 
random walk: 

dP( \ ] Q = k + P(i - 1; t) + k~P{% + 1; t) - {k + + k~) P(i; t), (13) 

(JjL 

where 

K^KoK-x 1 . . 



while 7 = k 3 and K 3 = k 3 /k 3 . 

From (13), the slow-dynamics stationary distribution obeys 



P*(i) = P* = constant. 

Observe that this distribution does not fulfill detailed balance since k + P*(i — 1) — k~P*(i) ^ 
0, unless k* = k~ . 

IV. THERMODYNAMIC STATE VARIABLES AND RELAXATION TO THE 
STATIONARY STATE 

Following [22H23], the enzyme internal energy, entropy, and Helmholtz free energy can 
be respectively defined as follows: 



u = -k B Tj2P(i,r,t)\ogP*(i,j), 

S = -k B J2 p (hr,t)^gP{i,j;t), 



F = U-TS = k B Tj2P(hj;t)\o$ 



P(i,r,t) 



P*(i,j) 

By substituting ^ into the above equations they can be rewritten as 

U = -A;BT^P(i;t)logP*(i)-A;flr2P(i;t)^P0'l*;*)log^(j|i): (15) 

i i j 

S = -A; i? ^P(z;t)logP(z;t)-fc B ^P(z;t)^P(jK;t)logP(j|z;t), (16) 

i i j 

F = ^T^P^^log^ + ^T^P^^^POI^^log^jl. (17) 

Observe that this way of writing the thermodynamic state variables renders a natural sep- 
aration of contributions from the fast and slow dynamics. 



The first term in the right hand side of Equation (17) is nothing else but the Kullback- 
Leibler divergence between distributions P{i;t) and P*(i) and so it is positive defined and 
only equals zero when the two distributions are identical. Similarly, the sum over j in the 
second term is the Kullback-Leibler divergence between P(j\i;t) and P*(j\i), it is positive 
defined, and only equals zero when P(j\i; t) = P*(j\i). From these considerations, the value 
of F can be used as an indicator of how far the system is from the stationary distribution. 



After differentiating Equation (17) and imposing the separation of time scales we obtain 

dF 
~dt 



Qkh — Tcr, 



(18) 



where Qhk is known as the housekeeping heat and is given by 



Q 



hk 



k B T J2 {Pfa ~ P{i + 1; t)k~) log — 



P*(i + l)k- 

2 

+ k B T ]T Pfo t) ( p 0'l*; t)K - p (j + t) k 7) log 

i j=l 



P*(j\i)k+ 
P*(j + l\i)k~ 



(19) 



while a is the entropy production rate: 



a = k B ^2(P(i;t)k + - P(i + 1; t)k~) log 



P{i-t)k+ 
P(i + l;t)k~ 



+ k B J2 Pfr t) E ~ p 0' + *)*7) i°s 



P(J\j;t)k+ 

P(j + l\i;t)k- 



(20) 



Observe that both and a have contributions from the slow (first term on the right hand 
side) and fast dynamics (second term). However, the fast dynamics contribution to Qhk 
vanishes because P*(j\i) complies with detailed balance, and so 



Qhk = k B T J2 {Pfr ' P(i + 1; t)k~) log — 



P*{i)k- 



(21) 



P*(i + l)k~ 

This result is in complete agreement with the interpretation of Qhk as the energy that has 
to be pumped into the system to drive the stationary state out from equilibrium (detailed 
balance) . 

In the present formalism, the enzyme molecule is implicitly assumed to be in equilibrium 
with a thermal bath. Concomitantly, the proper thermodynamic description is that of the 



Helmholtz free energy. As seen in Equation (17), F > and it only equals zero at the 
steady state. In other words, the value of F can be used as a measure of how distant the 
system state is from the stationary one, as we have previously asserted. Furthermore, it is 



not hard to prove from (18), (20), and (21) that dF/dt < 0, and that dF/dt = only when 
P(i]t) = P*(i) and P(j\i;t) = P*(J\i). Therefore, dF/dt can be interpreted as the rate of 
relaxation to the stationary distribution. 



We can see from (18), (20), and (21) that dF/dt can be decomposed into contributions 



from the slow and fast dymamics: 



dF 
~dt 



-^slow Pf&sti 
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(22) 



where 

f s1ow = k B J2 ^ + - P(i + 1; «)*-) log + (23) 

and 

2 P(?|ri)A; + 

F fast = * fl £ p(<; o £ (P01<; - ^(j + «)*7) log P(? - '' 1 , | /^ 7 - (24) 

i j=l ^ ' ' ' i 

The relaxation of the fast dynamics subspace to the corresponding quasi-stationary 
state — given by the slow-dynamics probability distribution — is determined by Pf as t- Accord- 
ing to the assumed separation of time scales, fast dynamics relaxation takes place without 
any noticeable change in the slow dynamics probability distribution (P(i;t)) and, in the 
long run, the slow dynamics are the ones that govern the system relaxation to the steady 
state. That is, the system relaxation rate is well approximated by P s iow 

Under the assumption that the probability distribution P(i; t) is slightly different from 
the stationary distribution P*\ P(i; t) = P* +e(i; t), with e(i; t) <C 1. Then, it is straightfor- 



ward to see from (23) that, in the neighborhood of the stationary distribution, the system 



relaxation rate (here denoted by £) is proportional to — see equation (14): 



K X K 2 K 3 -1 , . 



To better understand the result in Equation (25) let us analyze the significance of pa- 
rameters Ki (i = 1,2,3). From their definition, these parameters are nothing else but the 
association constants of the following reactions: Ep ^ E + P (Ki), E + S ^ E$ (K2), 
and Es ^ Ep (if 3), respectively. On the other hand, a reaction's association rate (K a) 
is related to its free energy change (AG) by Ka = exp(— AG/RT), with R the ideal gas 
constant and T the absolute temperature. These considerations allow us to visualize an 
enzymatic reaction as a process occurring along an energy profile like the one pictured in 
Figure [2} In such scheme, the enzyme states correspond to the local minima of the energy 
profile, and the transition probabilities are determined by the height of the energy barriers. 

The global free energy change (AGt) of an enzymatic reaction is given by AGt = AG% + 
AG2 + AG3 = —RTlog(KiK 2 K 3 ). Therefore, since the presence of an enzyme does not 
change AGt'- 

( AG T \ 
K 1 K 2 K 3 = exp I — -j^- J = constant. 

This restriction further implies that only two of the three Ki constants are independent. 
Without loss of generality we shall consider that K\ and K 2 are determined by the nature of 

8 



o 

Q- 

>> 

a> 

CD 
C 
LU 




— T" 

(M) 

Ep 



— I — 

(i,2) 
E+P 



i,3) 



0,4) 
E D 



Enzyme States 



FIG. 2: A cartoon representation of the Gibbs free energy profile for an enzymatic reaction. The 
minima correspond to the states the enzyme goes through, while the transition probabilities are 
determined by the height of the energy barriers. The individual free energy changes between 
adjacent states can be either positive or negative, but the global free energy change has to be 
negative in order for the reaction to proceed forward. 

the enzyme catalyzing the reaction, while K% = exp(—AG/RT)/KiK 2 . With this, Equation 

exp(-AG T /RT) - 1 



(25) can be rewritten as 



f oc 7- 



(26) 



1 + K x + K X K 2 

Note that the relaxation rate is a monotonic decreasing function of both K\ and K 2 . If we 
further take into consideration that Ki,K 2 > 0, it follows that the relaxation rate can be 
increased by making K\ and K 2 as close to zero as possible, and therefore by making AGi 
and AG 2 as positively large as possible. In particular, the maximum value of £ is attained 
when K\ = 0, regardless the value of K 2 . Thus, in order to increase the value of £ it is more 
important to increase the value of AGi than that of AG*2- 



On the other hand, we can see from Equation (26) that the relaxation rate is a monotonic 



decreasing function of AGt- That is, the more energetically favorable the global function 
is, the more stable the stationary distribution becomes. Similarly, £ is proportional to 
7 = — £3. Therefore, the system stationary distribution is more strongly stable when the 
substrate-to-product conversion process is more rapid. 
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V. REACTION VELOCITY AND EFFICIENCY OF ENERGY CONVERSION 



As we have seen, by exploiting the time-scale separation to simplify the reaction scheme, 
an enzyme can be modeled as a biased one-dimensional random walk with forward and 



backward transition probabilities k + and k~, respectively — see Equations (13) and (14). 
Therefore, by taking into consideration that the stationary probability distribution P* is 
constant, the reaction velocity can be defined as the forward minus the backward fluxes: 

*„4. , s *Gxp(-AG T /RT) - 1 . . 

» = ^-*-) = 7f T +Jfl + gl if, • (27) 

On the other hand, if we consider that a certain amount of energy is consumed during each 
forward step, and that this energy is waisted during backward steps, the system efficiency 
in the stationary state can be define as 

71 = 1 - k k^ = 1 -eM-^G T /RTy (28) 

Observe that v > is a decreasing function of AGt, and that v — when AGt = 0. In 
other words AGt < in order to have a positive reaction velocity. The efficiency r] is also a 
decreasing function of AGt, and rj = when AGt = 0. That is, both the reaction velocity 
and its efficiency can be increased by making AGt more negative. 



Note from equations (26) and (28) that rj oc £. Therefore, the discussion regarding the 
dependence of £ on K\ and K 2 applies as well to rj. In particular, we want to emphasize that 
the reaction velocity can be increased by making K\ and K 2 as close to zero as possible (and 
thus by making AGi and AG 2 as large as possible). However, varying K\ is more important 
since the maximum velocity can be achieved by setting K\ = 0, regardless the value of K 2 . 

Interestingly, the reaction efficiency is independent of K\ and K 2 . Thus, given that an 
enzyme does not alter the global free energy change of the reaction it catalyzes, this result 
implies that a reaction efficiency is the same regardless whether it is catalyzed or not. 



VI. CONCLUDING REMARKS 



In this work we have investigated the dynamic stability, as well as velocity and efficiency, 
of an enzymatic reaction with a rate limiting step, at the single molecule level. For this, 
we followed the ideas in [221 121], and used a recently developed formalism for performing 



10 



multiscale thermodynamic analysis on discrete-state, continuos-time, Markovian stochastic 
processes [22] • Our results can be summarized as follows: 



1. The dynamic and thermodynamic characteristics associated to the stationary proba- 
bility distribution are completely determined by the the Gibbs free energy changes of 
the enzymatic reaction steps: AGi (Ep ^ E + P), AG 2 (E + S ^ Es), and AG3 
(Es ^ Ep). Observe that the energies AGi (i = 1, 2, 3) are not all independent becase 
Y^=i AGi = AGt, and AGt (the global free energy change) is not modified by the 
enzyme. 

2. The stationary probability distribution is stable and the corresponding relaxation rate 
(£) is directly proportional to the global reaction velocity (v). 

3. Both £ and v are decreasing functions of AGt, and £, v — when AGt = 0. Thus, 
the global reaction accelerates and the stationary probability distribution turns more 
strongly stable as AGt is more negative. 

4. Both £ and v are increasing functions of AGi and AG 2 . The relaxation rate and the 
reaction velocity achieve their maximum value in the limit AGi — > 00, regardless the 
value of AG 2 . Contrarily, £ and v increase as AGi increases and converge to a value 
that depends on AGi as AG 2 — > 00. 

5. The efficiency (77) is a function of AGt, but it is independent of AGi and AG 2 - In 
particular, rj is a decreasing function of AGt, and 77 = when AGt = 0. That is, the 
reaction efficiency increases as AGt becomes more negative. 

Regarding the novelty of the above results, it is worth pointing out that, while some 
of them (like the one stating that v decreases as AGt increases) are well known, there are 
others (like those regarding the relaxation rate £) which are new to the best of our knowledge. 

On the other hand, when the separation of time scales is not possible, the system dynamic 
and thermodynamic characteristics will be determined by the rate constants kf and 
(i = 1,2,3), but the variables £, u, and rj, shall not necessarily be given in terms of the 



equilibrium dissociation constants Ki = k~ /kf, as in Equations (26)-(28). The reason being 
that the chemical reactions taking place between the synthesis of one product molecule and 
the next won't necessarily be close to chemical equilibrium. Some of the formerly enumerated 
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results may still be qualitatively true for enzymatic reactions without a rate limiting step, 
though. However, the study of such reactions is beyond the scope of the present work. In 
the forthcoming paragraphs we analyze the implications of the results listed above; in the 
understanding that, unless otherwise stated, these implications may not be exclusive for 
enzymatic reactions with a rate limiting step. 

Although in an enzymatic reaction the global free energy change is predetermined by the 
substrate, the products, and the nature of their reaction, the actual shape of the energy 
profile along the reaction coordinate depends on the enzyme structure. Our results confirm 
that although the efficiency is not affected by the shape of the energy profile, the reaction 
velocity and the strength of the stationary distribution stability can be highly improved by 
properly shaping this profile. This behavior is contrary to that observed in other systems 
(like thermal engines) in which variation in some parameters makes the velocity and the 
efficiency change in opposite directions. In other words, in this case no trade-off between 
efficiency and reaction velocity (and stability strength) need to be looked for, regarding 
the energy profile. We believe this may by one of the reasons why it has been possible 
for evolution to drive the structure of enzymes so the corresponding reaction velocity is 
increased by several orders of magnitude. 

Another feature worth noticing is the fact that the more energetically unfavorable reac- 
tions Ep f± E + P and E + S ^ E$ are, the faster the global enzymatic reaction is. This 
behavior can be understood by looking at Figure §2§. We can see there that large, positive 
AG i and AG2 values imply a very negative AG3. This further means that the backward 
reaction E$ ^— Ep is much less probable than the forward reaction E$ Ep. That is, 
the strategy to accelerate the global enzymatic reaction seems to be making the rate lim- 
iting step almost unidirectional, even though this implies that, in the rapid processes, the 
backward reactions are more probable than the corresponding forward reactions. Given the 
prominent role the rate limiting step takes in the above discussion, we believe it does not 
apply to other types of enzymatic reactions. 
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